####################################################################################################
# Replication Material for "How Political Careers affect Prime-Ministerial Performance" Table C1   #
# Authors: Florian Grotz, Ferdinand Müller-Rommel, Jan Berz, Corinna Kroeber, Marko Kukec          #
# Tested with R version 3.5.2                                                                      #
####################################################################################################

# This R code specifically replicates the calculation of individual cabinet and expert level ICC numbers
# It also contains an alternative estimation of ICC confidence intervals

library(foreign)
library(dplyr)

# please note that pre 0.18.0 versions of sjstats are required for the function get_var_re to calculate ICCs on the 
# individual expert and cabinet level
library(sjstats)

library(ICC)
library(readstata13)
library(arm)
library(lme4)
library(xtable)


# read the dataset
dat <- read.dta13("PMP_complete_expert_level.dta")

# save as R file
save(dat,file="PMP_dat.Rda")

load(file = "PMP_dat.Rda")

# check successful loading by number of cases
nrow(dat)


# two way random effects models
mconflict1 <- lmer(as.numeric(pmp_conflict1) ~ 1 +
                (1 | pmc_name) +
                (1 | expert), data = dat)

# icc function does not separate cabinet and expert level ICC, these will be calculated below
icc(mconflict1)

mconflict2 <- lmer(as.numeric(pmp_conflict2) ~ 1  + 
                     (1 | pmc_name) +
                     (1 | expert),
             data = dat)
icc(mconflict2)

mpolicy1 <- lmer(as.numeric(pmp_policy1) ~ 1  + 
                   (1 | pmc_name) +
                   (1 | expert), data = dat)
icc(mpolicy1)

mpolicy2 <- lmer(as.numeric(pmp_policy2) ~ 1  + 
                   (1 | pmc_name) +
                   (1 | expert), data = dat)
icc(mpolicy2)



mcrisis1 <- lmer(as.numeric(pmp_crisis1) ~ 1  + 
                   (1 | pmc_name) +
                   (1 | expert), data = dat)
icc(mcrisis1)


mcrisis2 <- lmer(as.numeric(pmp_crisis2) ~ 1  + 
                   (1 | pmc_name) +
                   (1 | expert), data = dat)
icc(mcrisis2)


mintern1 <- lmer(as.numeric(pmp_international1) ~ 1  + 
                   (1 | pmc_name) +
                   (1 | expert), data = dat)
icc(mintern1)

mintern2 <- lmer(as.numeric(pmp_international2) ~ 1  + 
                   (1 | pmc_name) +
                   (1 | expert), data = dat)
icc(mintern2)



mparl1 <- lmer(as.numeric(pmp_parliament1) ~ 1  + 
                 (1 | pmc_name) +
                 (1 | expert), data = dat)

icc(mparl1)

mparl2 <- lmer(as.numeric(pmp_parliament2) ~ 1  + 
                 (1 | pmc_name) +
                 (1 | expert), data = dat)

icc(mparl2)

mparty1 <- lmer(as.numeric(pmp_party1) ~ 1  + 
                  (1 | pmc_name) +
                  (1 | expert), data = dat)

icc(mparty1)

mparty2 <- lmer(as.numeric(pmp_party2) ~ 1  + 
                  (1 | pmc_name) +
                  (1 | expert), data = dat)

icc(mparty2)

mrating <- lmer(as.numeric(pmp_rating) ~ 1  + 
                  (1 | pmc_name) +
                  (1 | expert), data = dat)
icc(mrating)

# calculate the confidence intervals of the cabinet level ICC by bootstrapping
ci.mconflict1 <- bootMer(mconflict1, calc.icc, nsim=1000)
icc.ci.conflict1 <- quantile(ci.mconflict1$t, c(0.025, 0.975))

ci.mconflict2 <- bootMer(mconflict2, calc.icc, nsim=1000)
icc.ci.conflict2 <-quantile(ci.mconflict2$t, c(0.025, 0.975))

ci.mpolicy1 <- bootMer(mpolicy1, calc.icc, nsim=1000)
icc.ci.policy1 <-quantile(ci.mpolicy1$t, c(0.025, 0.975))

ci.mpolicy2 <- bootMer(mpolicy2, calc.icc, nsim=1000)
icc.ci.policy2 <-quantile(ci.mpolicy2$t, c(0.025, 0.975))

ci.mcrisis1 <- bootMer(mcrisis1, calc.icc, nsim=1000)
icc.ci.crisis1 <-quantile(ci.mcrisis1$t, c(0.025, 0.975))

ci.mcrisis2 <- bootMer(mcrisis2, calc.icc, nsim=1000)
icc.ci.crisis2 <-quantile(ci.mcrisis2$t, c(0.025, 0.975))

ci.mintern1 <- bootMer(mintern1, calc.icc, nsim=1000)
icc.ci.intern1 <-quantile(ci.mintern1$t, c(0.025, 0.975))

ci.mintern2 <- bootMer(mintern2, calc.icc, nsim=1000)
icc.ci.intern2 <-quantile(ci.mintern2$t, c(0.025, 0.975))

ci.mparl1 <- bootMer(mparl1, calc.icc, nsim=1000)
icc.ci.parl1 <-quantile(ci.mparl1$t, c(0.025, 0.975))

ci.mparl2 <- bootMer(mparl2, calc.icc, nsim=1000)
icc.ci.parl2 <-quantile(ci.mparl2$t, c(0.025, 0.975))

ci.mparty1 <- bootMer(mparty1, calc.icc, nsim=1000)
icc.ci.party1 <-quantile(ci.mparty1$t, c(0.025, 0.975))

ci.mparty2 <- bootMer(mparty2, calc.icc, nsim=1000)
icc.ci.party2 <-quantile(ci.mparty2$t, c(0.025, 0.975))

ci.mrating <- bootMer(mrating, calc.icc, nsim=1000)
icc.ci.rating <-quantile(ci.mrating$t, c(0.025, 0.975))

iccs_cabinet_cis <- rbind(icc.ci.conflict1,icc.ci.conflict2,icc.ci.policy1,icc.ci.policy2,
                          icc.ci.crisis1,icc.ci.crisis2,icc.ci.intern1,icc.ci.intern2,
                          icc.ci.parl1,icc.ci.parl2,icc.ci.party1,icc.ci.party2,icc.ci.rating)

# calculate specific cabinet levels iccs – separation requires out of icc() calculation
conflict1.icc <- get_re_var(mconflict1)[2] / (sum(get_re_var(mconflict1)) + get_re_var(mconflict1, "sigma_2")) 
conflict2.icc <- get_re_var(mconflict2)[2] / (sum(get_re_var(mconflict2)) + get_re_var(mconflict2, "sigma_2")) 
policy1.icc <- get_re_var(mpolicy1)[2] / (sum(get_re_var(mpolicy1)) + get_re_var(mpolicy1, "sigma_2")) 
policy2.icc <- get_re_var(mpolicy2)[2] / (sum(get_re_var(mpolicy2)) + get_re_var(mpolicy2, "sigma_2")) 
crisis1.icc <- get_re_var(mcrisis1)[2] / (sum(get_re_var(mcrisis1)) + get_re_var(mcrisis1, "sigma_2")) 
crisis2.icc <- get_re_var(mcrisis2)[2] / (sum(get_re_var(mcrisis2)) + get_re_var(mcrisis2, "sigma_2")) 
intern1.icc <- get_re_var(mintern1)[2] / (sum(get_re_var(mintern1)) + get_re_var(mintern1, "sigma_2")) 
intern2.icc <- get_re_var(mintern2)[2] / (sum(get_re_var(mintern2)) + get_re_var(mintern2, "sigma_2")) 
parl1.icc <- get_re_var(mparl1)[2] / (sum(get_re_var(mparl1)) + get_re_var(mparl1, "sigma_2")) 
parl2.icc <- get_re_var(mparl2)[2] / (sum(get_re_var(mparl2)) + get_re_var(mparl2, "sigma_2")) 
party1.icc <- get_re_var(mparty1)[2] / (sum(get_re_var(mparty1)) + get_re_var(mparty1, "sigma_2")) 
party2.icc <- get_re_var(mparty2)[2] / (sum(get_re_var(mparty2)) + get_re_var(mparty2, "sigma_2")) 
rating.icc <- get_re_var(mrating)[2] / (sum(get_re_var(mrating)) + get_re_var(mrating, "sigma_2")) 

iccs <- rbind(conflict1.icc,conflict2.icc,policy1.icc,policy2.icc,crisis1.icc,crisis2.icc,
              intern1.icc,intern2.icc,parl1.icc,parl2.icc,party1.icc,party2.icc, rating.icc)

# expert level specific iccs
conflict1.icc2 <- get_re_var(mconflict1)[1] / (sum(get_re_var(mconflict1)) + get_re_var(mconflict1, "sigma_2")) 
conflict2.icc2 <- get_re_var(mconflict2)[1] / (sum(get_re_var(mconflict2)) + get_re_var(mconflict2, "sigma_2")) 
policy1.icc2 <- get_re_var(mpolicy1)[1] / (sum(get_re_var(mpolicy1)) + get_re_var(mpolicy1, "sigma_2")) 
policy2.icc2 <- get_re_var(mpolicy2)[1] / (sum(get_re_var(mpolicy2)) + get_re_var(mpolicy2, "sigma_2")) 
crisis1.icc2 <- get_re_var(mcrisis1)[1] / (sum(get_re_var(mcrisis1)) + get_re_var(mcrisis1, "sigma_2")) 
crisis2.icc2 <- get_re_var(mcrisis2)[1] / (sum(get_re_var(mcrisis2)) + get_re_var(mcrisis2, "sigma_2")) 
intern1.icc2 <- get_re_var(mintern1)[1] / (sum(get_re_var(mintern1)) + get_re_var(mintern1, "sigma_2")) 
intern2.icc2 <- get_re_var(mintern2)[1] / (sum(get_re_var(mintern2)) + get_re_var(mintern2, "sigma_2")) 
parl1.icc2 <- get_re_var(mparl1)[1] / (sum(get_re_var(mparl1)) + get_re_var(mparl1, "sigma_2")) 
parl2.icc2 <- get_re_var(mparl2)[1] / (sum(get_re_var(mparl2)) + get_re_var(mparl2, "sigma_2")) 
party1.icc2 <- get_re_var(mparty1)[1] / (sum(get_re_var(mparty1)) + get_re_var(mparty1, "sigma_2")) 
party2.icc2 <- get_re_var(mparty2)[1] / (sum(get_re_var(mparty2)) + get_re_var(mparty2, "sigma_2")) 
rating.icc2 <- get_re_var(mrating)[1] / (sum(get_re_var(mrating)) + get_re_var(mrating, "sigma_2")) 

iccs2 <- rbind(conflict1.icc2,conflict2.icc2,policy1.icc2,policy2.icc2,crisis1.icc2,crisis2.icc2,
              intern1.icc2,intern2.icc2,parl1.icc2,parl2.icc2,party1.icc2,party2.icc2,rating.icc2)

allicc <- cbind(iccs,iccs2,iccs_cabinet_cis)

allicc <- round(allicc, 2)

write.table(allicc, file = "icc.txt", sep = ";", quote = FALSE, row.names = T)
